ggplot(Data, aes(x=Distance_um, y = Sr88_ppm, color = Habitat3)) + geom_point() +
facet_wrap(~Size*agency_id)
# Some fish are obviously anadromous and resident. Some are confusing - perhaps spending time in fjord?
ggplot(AData, aes(x=Distance_um, y = Sr88_ppm, color = Habitat3)) + geom_point() +
facet_wrap(~agency_id)
# fish not sure what phenotype they are as Sr increase but not to same concs as more 'obvious' sea trout with repeat migrations
Uncertain = c("LT21_3187", "LT_3302", "LT21_3241", "LT21_3302", "ST21_2006")
#Mark if adult fish ever went to sea in otolith material outside the core
Sea = subset(AData, Habitat3 =="Sea" & Distance_um > 250)
range(na.omit(unique(Sea$Age)))
Sea = unique(Sea$agency_id)
AData = AData %>% mutate(Phenotype = case_when(agency_id %in% Uncertain ~ "Uncertain",
agency_id %in% c(Sea) ~ "Anadromous",
TRUE ~ "Resident"),
# BASED ON LENS DATA  PUT ALL UNUSURE FISH AS ANADROMOUS
Phenotype2 = case_when(agency_id %in% c(Sea) ~ "Anadromous",
agency_id %in% Uncertain ~ "Anadromous",
TRUE ~ "Resident"),
# BASED PURELY ON SR CUTOFF LT21_3187 IS PUT AS RESIDENT
Phenotype_cutoff_based = case_when(agency_id %in% c(Sea) ~ "Anadromous",
TRUE ~ "Resident"))
phenotypes = distinct(AData, agency_id, Phenotype, Phenotype2, Phenotype_cutoff_based)
#Does Ba help to interpret those 4? Nope - even some of the residents have really low and consistent Ba
ggplot(AData, aes(x=Distance_um, y = Ba138_ppm, color = Habitat3)) + geom_point() +
facet_wrap(~Phenotype*agency_id)+ylim(0,8)
# age at first outmigration
first_marine = filter(AData, Phenotype2 == "Anadromous" & Habitat == "Sea" & Distance_um > 250) %>%
group_by(agency_id) %>% slice(which.min(Distance_um)) %>% mutate(first_sea_entry = "y")
AData = left_join(AData, first_marine)
#First marine migration and ages at sea
# Fig. 6 -----------------------------------------------------------------
ggplot(filter(AData, !is.na(Habitat3), !is.na(Year), Year<=8), aes(x=Year, y = `Sr.Ca`, color = Habitat3)) + geom_point() +
geom_hline(yintercept = 2.42, linetype = "dotted", color="turquoise4")+
geom_point(data = first_marine, aes(x=Year, y = Sr.Ca), color = "black", size = 3)+
geom_line(data = filter(AData, !is.na(Habitat3), !is.na(Year)), aes(x=Year, y = Sr.Ca), color = "black", size = 0.5, alpha=0.5)+
facet_wrap(~Phenotype*agency_id)+labs(x = "Estimated age", y = "Otolith Sr:Ca (mmol/mol)", color = "Assigned\nhabitat")+
scale_x_continuous(breaks = c(0:9))+
#xlim(0,8)+
theme_bw(base_size=15)+
theme(panel.grid.major.x = element_line(color = "grey80"),   # Style for gridlines
panel.grid.minor.x = element_blank(),
panel.grid.minor.y = element_blank(),
panel.grid.major.y = element_blank())+Theme +
scale_color_manual(values=c("#999999","#56B4E9"))                # Remove minor gridlines
ggsave('Plots/Fig. 6.eps', width = 12, height = 10, units = "in", dpi = 600)
mean(first_marine$Year[first_marine$Phenotype=="Anadromous"])
sd(first_marine$Year[first_marine$Phenotype=="Anadromous"])
# just look at annulus widths
inc_widths = filter(AData, !is.na(Increment) & Increment > 0)
ggplot(inc_widths, aes(x=Year, y = Increment, color = Phenotype)) +
geom_smooth(se=F) + geom_point() + geom_line(data = inc_widths, aes(group = agency_id), alpha = .3)
library(mgcv)  # ensure mgcv is loaded for gam support
growth <- ggplot(filter(inc_widths, !is.na(Habitat), Year <= 7),
aes(x = Year, y = Increment, color = Phenotype2)) +
geom_smooth(method = "gam", formula = y ~ s(x, bs = "cs",k=7), se = TRUE) +  # GAM smoother
geom_point(size = 4, alpha = 0.5, aes(shape = Habitat)) +
geom_line(data = filter(inc_widths, !is.na(Habitat), Year <= 7),
aes(group = agency_id), alpha = 0.3) +
labs(x = "Age", y = "Annulus width (µm)", color = "Phenotype") +
scale_x_continuous(breaks = c(1:7)) +
theme_bw(base_size = 12) +
theme(legend.position = "top") +
guides(color = guide_legend(nrow = 2),
shape = guide_legend(nrow = 2))
growth
# first try includes Uncertain fish but assumes ALL are anadromous (i.e. goes against threshold for LT21_3187)
# exclude rows where Year is not a whole number or for ages > 7 as only 1 fish
inc_widths2 = filter(inc_widths, !is.na(Habitat), Year<=7) %>%
mutate(Phenotype2 = as.factor(Phenotype2),
Phenotype_cutoff_based = as.factor(Phenotype_cutoff_based),
agency_id = as.factor(agency_id))
mod = gam(Increment ~ Phenotype2 + Habitat + # is there an overall difference in growth among phenotypes or habitats - Nope
#s(Year, k=6)+ # Exclude as if this is included all the variation is explained here (ie ontogeny)
s(Year, by=Phenotype2, k=6)+# is there a different slope between phenotypes? Some weak suggestion of yes - for age 1 only
s(agency_id, Year, bs = 're'), # random effect of animal
data=inc_widths2, na.action = "na.omit", method="REML")
summary(mod)
library(gratia)
diff1 <- difference_smooths(mod, smooth = "s(Year)")
## Fig 7b ------------------------------------------------------------------
gam_diff = draw(diff1) & theme_bw() & xlab("Age") & scale_x_continuous(breaks = c(1:7)); gam_diff
cowplot::plot_grid(growth, gam_diff, ncol=1, align = "v", labels = c("a","b"))
ggsave("Plots/Fig. 7.eps", height = 7, width=5,dpi=600)
##Remove strange fish
# Try and excludes Uncertain fish - sensitivity analysis
inc_widths3 = filter(inc_widths2, !Phenotype=="Uncertain")
mod2 = gam(Increment ~ Phenotype2 + Habitat + # is there an overall difference in growth among phenotypes or habitats - Nope
#s(Year, k=6)+ # Exclude as if this is included all the variation is explained here (ie ontogeny)
s(Year, by=Phenotype2, k=6)+# is there a different slope between phenotypes? Some weak suggestion of yes - for age 1 only
s(agency_id, Year, bs = 're'), # random effect of animal
data=inc_widths3, na.action = "na.omit", method="REML")
summary(mod2)
diff2 <- difference_smooths(mod2, smooth = "s(Year)")
# #Fig. S8 -----------------------------------------------------------------
gam_diff2 = draw(diff2) & theme_bw()&theme(panel.grid.minor = element_blank(),axis.title=element_text(size=20),axis.text=element_text(size=18)) & xlab("Age") & scale_x_continuous(breaks = c(1:7)); gam_diff2
gam_diff2
# Difference between smooths is sig for age 1 only (i.e. Anadromous fish have faster first year growth but caveat that we have 3 resident fish!!)
ggsave("Plots/Fig. S5.eps", height = 5, width=5,dpi=600)
EData<- read_csv("Data Outputs/alleyedata.csv")
EData=merge(EData, phenotypes, by=c("agency_id"))
EData <- EData %>%
mutate(
LineTypeGroup = ifelse(agency_id == "LT21_3187", "LT21_3187", "other")
)
ggplot(EData, aes(x = diameter_um, y = d34S,
color = Phenotype,
group = agency_id,
linetype = LineTypeGroup)) +
facet_wrap(~Phenotype) +
geom_line(size = 1) +
theme_bw() + Theme +
xlab("Eye lens diameter") +
ylab(expression(paste("Lens ", delta^34, "S"))) +
theme(strip.text = element_text(size = 16)) +
scale_linetype_manual(values = c("LT21_3187" = "dashed", "other" = "solid"))
## Fig. S7 d34S ------------------------------------------------------------
ggsave("Plots/Fig. S4.eps", height = 6, width=9,units = "in",dpi=600)
#General Data Viz
#Bring in required packages
library(readxl)
library(tidyr)
library(plyr)
library(dplyr)
library(ggplot2)
library(oce)
library(segmented)
library(zoo)
library(mixtools)
library(randomForest)
library(here)
library(caret)
library(tidyverse)
# Load Packages -----------------------------------------------------------
#Theme
Theme = theme(panel.grid.major = element_blank(),panel.grid.minor = element_blank(),axis.title=element_text(size=20),axis.text=element_text(size=18))
# Data Prep ----------------------------------------------------
EData = read.csv("Data Outputs/processedotodata2_DS.csv")
#Subset Juvaniles
JD=subset(EData, Size=="Juv")
range(na.omit(unique(JD$LRadius)))
mean(na.omit(unique(JD$LRadius)))
#Subset Adults
AData = subset(EData, Size =="AD")
# Example Tracer profiles -------------------------------------------------
####Select example Sr,Ba, Zn plot
Ex = subset(AData, agency_id =="LT21_3196")
Ex= mutate(Ex, `Ba.Ca` = (Ba138_ppm/137.327)/((Ca43_ppm/40.078)/1000000))
Ex= mutate(Ex, `Mg.Ca` = (Mg24_ppm/24.305)/((Ca43_ppm/40.078)/1000000))
Ex= mutate(Ex, `Zn.Ca` = (Zn66_ppm/65.38)/((Ca43_ppm/40.078)/1000000))
Sr = ggplot(Ex, aes(x=Distance_um/1000, y=`Sr.Ca`))+geom_line(size=0.8)+theme_bw()+Theme+ylab("Sr:Ca (mmol/mol)")+
theme(axis.title.x=element_blank(),
axis.text.x=element_blank(),
axis.ticks.x=element_blank())+ theme(plot.margin = unit(c(0.1, 0.1, 0.1, 0.1), "cm"))#+ xlab("Estimated Length (mm)")
Ba= ggplot(Ex, aes(x=Distance_um/1000, y=`Ba.Ca`))+geom_line(size=0.8)+theme_bw()+Theme+ylab(expression(paste("Ba:Ca (",mu,"mol/mol)")))+
theme(axis.title.x=element_blank(),
axis.text.x=element_blank(),
axis.ticks.x=element_blank())+ theme(plot.margin = unit(c(0.1, 0.1, 0.1, 0.1), "cm"))#+ xlab("Estimated Length (mm)")
Mg = ggplot(Ex, aes(x=Distance_um/1000, y=`Mg.Ca`))+geom_line(size=0.8)+theme_bw()+Theme+ xlab(expression(paste("Otolith Radius (mm)")))+ylab(expression(paste("Mg:Ca (",mu,"mol/mol)")))+ theme(plot.margin = unit(c(0.1, 0.1, 0.1, 0.1), "cm"))
Zn = ggplot((Ex)[!is.na(Ex$`Zn.Ca`),], aes(x=(Distance_um/1000), y=(`Zn.Ca`)))+geom_line(size=0.8)+theme_bw()+Theme+ xlab(expression(paste("Otolith Radius (mm)")))+ylab(expression(paste("Zn:Ca (",mu,"mol/mol)")))+ theme(plot.margin = unit(c(0.1, 0.1, 0.1, 0.1), "cm"))
#Sulphur
EData =  read_csv("Data Outputs/alleyedata.csv")
LEeye = subset(EData, num ==1)
ggplot(LEeye, aes(x=diameter_um, y=length_mm))+geom_point()+geom_smooth(method="lm")
m1=lm(length_mm~diameter_um, LEeye)
summary(m1)
m1$coefficients[2]
Ex2 = subset(EData, agency_id =="LT21_3196")
#Adjust eye lens diameter to radius
#Bring in cores
Ex2=dplyr::select(Ex2,d13C,d34S,diameter_um,agency_id,num,tag)
#Adjust
Ex2$Rad=Ex2$diameter_um/2
Ex2 = Ex2%>%
mutate(Rad2 = (Rad+lead(Rad))/2)
Sl = ggplot(Ex2, aes(x=Rad2/1000, y=d34S))+geom_line(size=0.8)+theme_bw()+Theme+
xlab(expression(paste("Lens Radius (mm)")))+geom_point()+ylim(5,11)+xlim(0,1.8)+
scale_y_continuous( limits = c(4.5, 11),sec.axis = sec_axis(trans = ~ . -35,name = expression(paste("Lens ",delta^13, "C")) ))+geom_point(aes(y=d13C+35),color="red")+
geom_line(aes(y=d13C+35),size=0.8,color="red")+
ylab(expression(paste("Lens ",delta^34, "S")))+ theme(plot.margin = unit(c(0.1, 0.1, 0.1, 0.1), "cm"))
blank=ggplot(Ex, aes(x=Distance_um/1000, y=))+theme_bw()+Theme+ xlab(expression(paste("Otolith Radius (mm)")))+ theme(plot.margin = unit(c(0.1, 0.1, 0.1, 0.1), "cm"))
## Fig. 3 ------------------------------------------------------------------
F3 =plot_grid(Sr, Ba, Sl,Zn, Mg,blank, ncol=3,labels=c("a","b","c","d","e","f"), align = "v")
## Fig. 3 ------------------------------------------------------------------
F3 =plot_grid(Sr, Ba, Sl,Zn, Mg,blank, ncol=3,labels=c("a","b","c","d","e","f"), align = "v")
F3
ggsave("Plots/Fig. 3.tiff", plot = F3, width = 15, height = 12, units = "in", dpi = 600)
# Compare Chemical vs Visual Ageing ----------------------------------------
#Chem age
Data = read.csv("Data Outputs/prepedotodata_DS.csv")
tag = data.frame("Sheet"=unique(Data$Sheet), "agency_id"=unique(Data$agency_id))
myFiles <- list.files("Raw Data/Chem Ages")
#Chem Age
results=NULL
for(i in levels(as.factor(unique(Data$agency_id)))){tryCatch({
AgData = read.csv(paste0("Raw Data/Chem Ages/",i))
AgData$agency_id = i
AgData$Year = NA
AgData$Age = nrow(AgData)
for(t in c(1:nrow(AgData))){
AgData[t,]$Year = t
}
results = rbind(results, AgData)}, error=function(e){})
}
results=merge(results, tag, by=c("agency_id"))
results=results%>% dplyr::select(Age, Sheet,agency_id)
results=unique(results)
# Compare Chemical vs Visual Ageing ----------------------------------------
#Import
Chem=results
Vis = read.csv("Data Outputs/processedOtoData_DS.csv")
Reread = read_excel("Data Outputs/Oto Rereads.xlsx")
#Cut to just age
Vis = Vis%>%
dplyr::select(Sheet, Age)%>%
na.omit()%>%
unique()
#Merge Reareads
Chem=merge(unique(Chem), Reread, by=c("Sheet"),all = TRUE)
Chem=subset(Chem,Reader =="Anna")
Vis=merge(unique(Vis), Reread, by=c("Sheet"),all = TRUE)
Vis=subset(Vis,Reader =="Coralie")
#See diferences
Chem$Diff=(Chem$Age.y-Chem$Age.x)
Chem$Type="Chemical"
Vis$Diff=(Vis$Age.y-Vis$Age.x)
Vis$Type="Visual"
#Merge
Comp = rbind(dplyr::select(Vis, Diff, Type),dplyr::select(Chem, Diff, Type))
#Plot
S2 = ggplot(Comp, aes(y=abs(Diff), x=Type))+geom_boxplot(outliers = FALSE)+geom_jitter(width=0.05,height=0.01,alpha=0.5,size=2)+ylim(-1,7)+theme_bw()+Theme+ylab("Difference in assigned age")+xlab("Ageing method")
ggsave("Plots/Fig. S2.eps", plot = S2, width = 5, height = 4, units = "in", dpi = 600)
## Fig. 3 ------------------------------------------------------------------
F3 =plot_grid(Sr, Ba, Sl,Zn, Mg,blank, ncol=3,labels=c("a","b","c","d","e","f"), align = "v")
F3
##Data analysis of growth between habitats
# Load libraries ----------------------------------------------------------
library(readxl)
library(tidyr)
library(plyr)
library(dplyr)
library(ggplot2)
library(oce)
library(segmented)
library(zoo)
library(mixtools)
library(randomForest)
library(here)
library(caret)
library(tidyverse)
library(mgcv)
#Theme
Theme = theme(panel.grid.major = element_blank(),panel.grid.minor = element_blank(),axis.title=element_text(size=20),axis.text=element_text(size=18))
# Import Data-------------------------------------------------------------------------
Data = read.csv("Data Outputs/processedotodata2_DS.csv")
AData = subset(Data, Size  =="AD")
JData = subset(Data, Size  =="Juv")
# Data Vis and Exploration ------------------------------------------------
# Threshold for marine migration based on max Sr value observed in immature reference samples outside of the core (250um otolith radius)
max(JData$Sr88_ppm[JData$Distance_um > 250], na.rm = T)
max(JData$Sr.Ca[JData$Distance_um > 250], na.rm = T)
ggplot(Data, aes(x=Distance_um, y = Sr88_ppm, color = Habitat3)) + geom_point() +
facet_wrap(~Size*agency_id)
# Some fish are obviously anadromous and resident. Some are confusing - perhaps spending time in fjord?
ggplot(AData, aes(x=Distance_um, y = Sr88_ppm, color = Habitat3)) + geom_point() +
facet_wrap(~agency_id)
# fish not sure what phenotype they are as Sr increase but not to same concs as more 'obvious' sea trout with repeat migrations
Uncertain = c("LT21_3187", "LT_3302", "LT21_3241", "LT21_3302", "ST21_2006")
#Mark if adult fish ever went to sea in otolith material outside the core
Sea = subset(AData, Habitat3 =="Sea" & Distance_um > 250)
range(na.omit(unique(Sea$Age)))
Sea = unique(Sea$agency_id)
AData = AData %>% mutate(Phenotype = case_when(agency_id %in% Uncertain ~ "Uncertain",
agency_id %in% c(Sea) ~ "Anadromous",
TRUE ~ "Resident"),
# BASED ON LENS DATA  PUT ALL UNUSURE FISH AS ANADROMOUS
Phenotype2 = case_when(agency_id %in% c(Sea) ~ "Anadromous",
agency_id %in% Uncertain ~ "Anadromous",
TRUE ~ "Resident"),
# BASED PURELY ON SR CUTOFF LT21_3187 IS PUT AS RESIDENT
Phenotype_cutoff_based = case_when(agency_id %in% c(Sea) ~ "Anadromous",
TRUE ~ "Resident"))
phenotypes = distinct(AData, agency_id, Phenotype, Phenotype2, Phenotype_cutoff_based)
# age at first outmigration
first_marine = filter(AData, Phenotype2 == "Anadromous" & Habitat == "Sea" & Distance_um > 250) %>%
group_by(agency_id) %>% slice(which.min(Distance_um)) %>% mutate(first_sea_entry = "y")
AData = left_join(AData, first_marine)
ggplot(filter(AData, !is.na(Habitat3), !is.na(Year), Year<=8), aes(x=Year, y = `Sr.Ca`, color = Habitat3)) + geom_point() +
geom_hline(yintercept = 2.42, linetype = "dotted", color="turquoise4")+
geom_point(data = first_marine, aes(x=Year, y = Sr.Ca), color = "black", size = 3)+
geom_line(data = filter(AData, !is.na(Habitat3), !is.na(Year)), aes(x=Year, y = Sr.Ca), color = "black", size = 0.5, alpha=0.5)+
facet_wrap(~Phenotype*agency_id)+labs(x = "Estimated age", y = "Otolith Sr:Ca (mmol/mol)", color = "Assigned\nhabitat")+
scale_x_continuous(breaks = c(0:9))+
#xlim(0,8)+
theme_bw(base_size=15)+
theme(panel.grid.major.x = element_line(color = "grey80"),   # Style for gridlines
panel.grid.minor.x = element_blank(),
panel.grid.minor.y = element_blank(),
panel.grid.major.y = element_blank())+Theme +
scale_color_manual(values=c("#999999","#56B4E9"))                # Remove minor gridlines
ggplot(filter(AData, !is.na(Habitat3), !is.na(Year), Year<=8), aes(x=Year, y = `Sr.Ca`, color = Habitat3)) + geom_point() +
geom_hline(yintercept = 2.42, linetype = "dotted", color="turquoise4")+
geom_point(data = first_marine, aes(x=Year, y = Sr.Ca), color = "black", size = 3)+
geom_line(data = filter(AData, !is.na(Habitat3), !is.na(Year)), aes(x=Year, y = Sr.Ca), color = "black", size = 0.5)+
facet_wrap(~Phenotype*agency_id)+labs(x = "Estimated age", y = "Otolith Sr:Ca (mmol/mol)", color = "Assigned\nhabitat")+
scale_x_continuous(breaks = c(0:9))+
#xlim(0,8)+
theme_bw(base_size=15)+
theme(panel.grid.major.x = element_line(color = "grey80"),   # Style for gridlines
panel.grid.minor.x = element_blank(),
panel.grid.minor.y = element_blank(),
panel.grid.major.y = element_blank())+Theme +
scale_color_manual(values=c("#999999","#56B4E9"))                # Remove minor gridlines
ggsave('Plots/Fig. 6.eps', width = 12, height = 10, units = "in", dpi = 600)
mean(first_marine$Year[first_marine$Phenotype=="Anadromous"])
sd(first_marine$Year[first_marine$Phenotype=="Anadromous"])
# just look at annulus widths
inc_widths = filter(AData, !is.na(Increment) & Increment > 0)
ggplot(inc_widths, aes(x=Year, y = Increment, color = Phenotype)) +
geom_smooth(se=F) + geom_point() + geom_line(data = inc_widths, aes(group = agency_id), alpha = .3)
library(mgcv)  # ensure mgcv is loaded for gam support
growth <- ggplot(filter(inc_widths, !is.na(Habitat), Year <= 7),
aes(x = Year, y = Increment, color = Phenotype2)) +
geom_smooth(method = "gam", formula = y ~ s(x, bs = "cs",k=7), se = TRUE) +  # GAM smoother
geom_point(size = 4, alpha = 0.5, aes(shape = Habitat)) +
geom_line(data = filter(inc_widths, !is.na(Habitat), Year <= 7),
aes(group = agency_id), alpha = 0.3) +
labs(x = "Age", y = "Annulus width (µm)", color = "Phenotype") +
scale_x_continuous(breaks = c(1:7)) +
theme_bw(base_size = 12) +
theme(legend.position = "top") +
guides(color = guide_legend(nrow = 2),
shape = guide_legend(nrow = 2))
growth
# first try includes Uncertain fish but assumes ALL are anadromous (i.e. goes against threshold for LT21_3187)
# exclude rows where Year is not a whole number or for ages > 7 as only 1 fish
inc_widths2 = filter(inc_widths, !is.na(Habitat), Year<=7) %>%
mutate(Phenotype2 = as.factor(Phenotype2),
Phenotype_cutoff_based = as.factor(Phenotype_cutoff_based),
agency_id = as.factor(agency_id))
mod = gam(Increment ~ Phenotype2 + Habitat + # is there an overall difference in growth among phenotypes or habitats - Nope
#s(Year, k=6)+ # Exclude as if this is included all the variation is explained here (ie ontogeny)
s(Year, by=Phenotype2, k=6)+# is there a different slope between phenotypes? Some weak suggestion of yes - for age 1 only
s(agency_id, Year, bs = 're'), # random effect of animal
data=inc_widths2, na.action = "na.omit", method="REML")
summary(mod)
library(gratia)
diff1 <- difference_smooths(mod, smooth = "s(Year)")
## Fig 7b ------------------------------------------------------------------
gam_diff = draw(diff1) & theme_bw() & xlab("Age") & scale_x_continuous(breaks = c(1:7)); gam_diff
cowplot::plot_grid(growth, gam_diff, ncol=1, align = "v", labels = c("a","b"))
ggsave("Plots/Fig. 7.tiff", height = 7, width=5,dpi=600)
#Freshwater Movement Comparison
# Load Packages -----------------------------------------------------------
library(readxl)
library(tidyr)
library(plyr)
library(dplyr)
library(ggplot2)
library(oce)
library(segmented)
library(zoo)
library(mixtools)
library(randomForest)
library(here)
library(caret)
library(tidyverse)
#Theme
Theme = theme(panel.grid.major = element_blank(),panel.grid.minor = element_blank(),axis.title=element_text(size=20),axis.text=element_text(size=18))
# Visualise Eye Data ------------------------------------------------------
EData = read.csv("Data Outputs/alleyedata.csv")
str(EData)
#See what period sampling
Samp <- EData %>%
group_by(agency_id) %>%
arrange(num, .by_group = TRUE) %>%   # use the variable that defines order
mutate(diameter_diff_um = diameter_um - lag(diameter_um)) %>%
ungroup()
Samp <- Samp %>%
filter(grepl("LAYER", tag)) %>%
group_by(agency_id) %>%
arrange(num, .by_group = TRUE) %>%   # change if a different order defines “first”
slice(2) %>%
ungroup()
Samp=subset(Samp,Size=="Juv")
mean(na.omit(Samp$diameter_diff_um))/2 #to get radius dif
2*sd(na.omit(Samp$diameter_diff_um))/2
#Core size
Core = subset(EData, tag=="CORE")
#Subset to most recent growth and remove non paired samples
Lakes = subset(EData, tag =="JELLY" &Lake != "HEGA" & Size == "Juv"&agency_id!="LT21_3074"&agency_id!="ST21_2055"&agency_id!="ST21_2062"&agency_id!="ST21_2260"&agency_id!="ST21_2280")
range(na.omit(Lakes$diameter_um))
Lakes = Lakes[,c(6,7,8,12,24,22)]
str(Lakes)
library(data.table)
LakesL1 <- melt(setDT(Lakes), id.vars = c("agency_id", "Lake", "Size"), variable.name = "Element")
#Visualise differences in stable isootpes between lakes
Lakeslt = subset(Lakes, Lake == "LT")
Lakesst = subset(Lakes, Lake == "ST")
Lakes = mutate(Lakes, Lake = recode(Lake,ST = "Storvatnet", LT = "Litjvatnet"))
#Individual Plots
B1 = ggplot(Lakes, aes(x=Lake, y=d13C))+geom_boxplot(outliers = FALSE)+geom_jitter(width=0.04,size=2,alpha=0.5)+theme_bw()+Theme+xlab("Lake of capture")+
ylab(expression(paste("Lens ",delta^13, "C")))
wilcox.test(Lakeslt$d13C, Lakesst$d13C)
B2 = ggplot(Lakes, aes(x=Lake, y=d15N))+geom_boxplot(outliers = FALSE)+geom_jitter(width=0.04,size=2,alpha=0.5)+theme_bw()+Theme+xlab("")+
ylab(expression(paste("Lens ",delta^15, "N")))
wilcox.test(Lakeslt$d15N, Lakesst$d15N)
B3 = ggplot(Lakes, aes(x=Lake, y=d34S))+geom_boxplot(outliers = FALSE)+geom_jitter(width=0.04,size=2,alpha=0.5)+theme_bw()+Theme+xlab("")+
ylab(expression(paste("Lens ",delta^34, "S")))
wilcox.test(Lakeslt$d34S, Lakesst$d34S)
library(cowplot)
#Vis
plot_grid(B2,B1,B3,nrow=1,labels=c("(a)","(b)","(c)"))
# Visualise Otoliths ------------------------------------------------------
OData = read.csv("Data Outputs/processedotodata2_DS.csv")
OLakes= subset(OData, Lake != "HEGA")
#Check profiles
df = OLakes %>% pivot_longer(cols = c(14:22), names_to = "Element", values_to = "Conc_ppm") %>%
filter(!Element=="Ca43_ppm")
ggplot(filter(df), aes(Distance_um, Conc_ppm, group=agency_id)) +
geom_line(alpha=.5, size=.8)+
facet_wrap(~Element, scales = "free")
#Subset to juveniles
OLakes= subset(OData, Lake != "HEGA"&Size=="Juv"&agency_id!="ST21_2230")
#Subset to last 10 values for edge signal
OLakes2=OLakes%>%
group_by(agency_id)%>%
slice_tail(n=10)
colnames(OLakes2)
#Select one row with a mean value of these 10
OLakes2 = OLakes2%>%
dplyr::select(c(agency_id,Size,Lake,Na23_ppm,Mg24_ppm,P31_ppm,K39_ppm,Mn55_ppm,Zn66_ppm,Sr88_ppm,Ba138_ppm))#Remove calcium as fixed
OLakes2 <- OLakes2 %>%
group_by(agency_id) %>%
summarise(Size=first(Size), Lake=first(Lake),across(c(Na23_ppm, Mg24_ppm, P31_ppm,
K39_ppm, Mn55_ppm, Zn66_ppm,
Sr88_ppm, Ba138_ppm),
\(x) mean(x, na.rm = TRUE)))
OLakes2%>%
group_by(Lake)%>%
summarise(median(Mn55_ppm))
#Log
cols_to_log <- c("Na23_ppm", "Mg24_ppm", "P31_ppm", "K39_ppm",
"Mn55_ppm", "Zn66_ppm", "Sr88_ppm", "Ba138_ppm")
for (col in cols_to_log) {
OLakes2[[paste0(col, "_log")]] <- log10(OLakes2[[col]])
}
OLakes2 <- OLakes2[ , !(names(OLakes2) %in% cols_to_log)]
#Turn into long format
library(data.table)
LakesL <- melt(setDT(OLakes2), id.vars = c("agency_id", "Lake", "Size"), variable.name = "Element")
LakesL=subset(LakesL, Element != ("Li7")& Element !=  ("Ca43")&Element!= ("P31_ppm_log"))
LakesL = mutate(LakesL, Lake = recode(Lake, ST = "Storvatnet", LT = "Litjvatnet"))
LakesL$Element <- gsub("_ppm_log$", "", LakesL$Element)
LakesL$Element <- sub(
"^([A-Za-z]+)([0-9]+)$",
"\\2\\1",
LakesL$Element
)
#Vis
ggplot(LakesL, aes(x=Lake, y=(value)))+geom_boxplot(outliers = FALSE)+geom_jitter(width=0.04,size=2,alpha=0.5)+theme_classic()+Theme+
facet_wrap(~Element, scales = "free",labeller = labeller(Element = function(x) gsub("_ppm_log", "", x)))+
xlab("Lake of capture")+ylab("log Otolith trace\nelement conentration (ppm)")+theme(strip.text = element_text(size = 16))
#ggsave("Plots/Fig. S4.tiff", plot = S4, width = 16, height = 12, units = "in", dpi = 100)
# Combined ----------------------------------------------------------------
LakesL$Element <- factor(LakesL$Element, levels = unique(LakesL$Element))
LakesL=rbind(LakesL, LakesL1)
lab_fun <- function(x) {
idx <- match(x, levels(LakesL$Element))
# remove leading isotope numbers (e.g. 55Mn -> Mn)
x_no_num <- gsub("^[0-9]+", "", x)
# convert d15N -> delta^15*"N"
x_delta <- gsub("^d([0-9]+)(.*)", "delta^\\1*\"\\2\"", x_no_num)
ifelse(
idx <= 7,
paste0(
"Ln~(\"Otolith\"~\"", x_no_num, ":Ca\")"
),
paste0(
"\"Eye lens\"~", x_delta
)
)
}
ggplot(LakesL, aes(x = Lake, y = value)) +
geom_boxplot(outliers = FALSE) +
geom_jitter(width = 0.04, size = 2, alpha = 0.5) +
theme_classic() +
Theme +
facet_wrap(
~Element,
scales = "free",
ncol = 5,
labeller = as_labeller(lab_fun, default = label_parsed)
) +
xlab("Lake of capture") +
ylab("Tracer value") +
theme(strip.text = element_text(size = 16))
ggplot(LakesL, aes(x = Lake, y = value)) +
geom_boxplot(outliers = FALSE) +
geom_jitter(width = 0.04, size = 2) +
theme_classic() +
Theme +
facet_wrap(
~Element,
scales = "free",
ncol = 5,
labeller = as_labeller(lab_fun, default = label_parsed)
) +
xlab("Lake of capture") +
ylab("Tracer value") +
theme(strip.text = element_text(size = 16))
#Save
ggsave("Plots/Fig. 4.eps", width = 16, height = 10, units = "in", dpi = 600)
